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We provide a novel setup for generalizing the two-dimensional pseudospin 5 = 1/2 Dirac equation, 
arising in graphene's honeycomb lattice, to general pseudospin-S*. We engineer these band structures 
as a nearest-neighbor hopping Hamiltonian involving stacked triangular lattices. We obtain multi- 
layered low energy excitations around half-filling described by a two-dimensional Dirac equation of 
the form H = vpS ■ p, where S represents an arbitrary spin-S* (integer or half-integer). For integer- 
s', a fiat band appears, whose presence modifies qualitatively the response of the system. Among 
physical observables, the density of states, the optical conductivity and the peculiarities of Klein 
tunneling are investigated. We also study Chern numbers as well as the zero-energy Landau level 
degeneracy. By changing the stacking pattern, the topological properties are altered significantly, 
with no obvious analogue in multilayer graphene stacks. 
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I. INTRODUCTION 

Since the seminal work on monolayer graphene, a sin- 
gle sheet of carbon atoms forming a honeycomb lat- 
tice, in 2004i, a lot of attention has been focused on 
this material. Its low energy properties close to half 
filling (i.e. pristine graphene) are well described by 
a two-dimensional massless Dirac equation, with the 
speed of light replaced by the appropriate Fermi veloc- 
ity ^ 10^ m/s. Most of the unusual electronic prop- 
erties of this material can be traced back to the mass- 
less Dirac nature of its quasiparticles and their unusual 
Berry phase. These include its linearly vanishing density 
of states (DOS) around half filling resembling a d-wave 
superconductor, unusual Landau quantization in a per- 
pendicular magnetic field and the anomalous half-integer 
quantum Hall effect. Additionally, phenomena such as 
the universal optical conductivity and high optical trans- 
parency, Klein tunneling through electric barriers arc also 
distinguishing features. 

The appearance of the massless Dirac equation has 
triggered further research to find out whether other sys- 
tems can possess similar behaviour or even generaliza- 
tions of the 5 = 1/2 Dirac physics to e.g. higher di- 
mensions, including additional terms. In the context of 
ultracold atom in optical lattices, several proposals have 
been put forward to realize a generalization of graphene 
physics in terms of the S = 1 Dirac equation^Ti^. Gener- 
alizations to higher S with spin-dependent hoppings^ as 
well as with artificial magnetic field^ are also possible. 

Here, we present a family of lattices whose low energy 
excitations around given fillings are described by a gen- 
eralized two-dimensional Dirac equation, 



H = vpS ■ p. 



(1) 



where p — {px,Py,0), and S = {Sx, Sy, Sz) is the ma- 
trix representation of an arbitrary spin 5* (integer or half 



integer), and vp is the Fermi velocity. These lattices con- 
sist of stackings of triangular layers, and include slabs of 
face-centred cubic and hexagonal close-packed lattices as 
special cases. Technically, the notion Weyl Hamiltonian^ 
is more appropriate for Eq. ([1]) for S > 1/2, though 
wc refer to it as generalized Dirac equation (sometimes 
omitting the "generalized") since our motivation comes 
primarily from graphene and its pseudospin- 1/2 Dirac 
equation. 

Our model is characterized by considerable simplicity 
and tunability. Furthermore, it contains a unique feature, 
absent from previous lattice realizations of higher spin S 
Hamiltonians: the possibility of — through a simple lat- 
eral shift in the layer positions — changing the chiral prop- 
erties of individual intcrlaycr hoppings, without changing 
the spectrum. However, there arc considerable changes 
on other properties, for example topological properties 
such as the multiplicity of the zero-energy Landau level 
degeneracy in magnetic field. 

The paper is organized as follows. In Sec. II wc intro- 
duce the lattice and discuss some of its general features. 
We then proceed to analyze the properties resulting from 
such a band structure: density of states (Sec. HI.), op- 
tical conductivity (Sec. IV.), Chern numbers and spin 
Chcrn numbers of the band structure (Sec. V.), zero 
mode degeneracy in both a uniform as well as a nonuni- 
form magnetic field (Sec. VI.) and Klein tunneling (Sec. 
VII.) for the spin-1 case, focusing on tunneling into the 
flat band. We also derive the general matching condition 
for the wavefunction for arbitrary pseudospin-S'. Finally, 
the relevant symmetry properties are highlighted in an 
appendix. 

The 5 = 1/2 version is realized in graphene^ and on the 
surface of 3D topological insulators^. Of the recent pro- 
posals for the 5 = 1 case^"— , our construction includes the 
dice lattice. Our work in many respect is complementary 
to Ref. 0, where diverse properties (topology, transport) 
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of the spin-S* Dirac equation were studied using a differ- 
ent lattice realization with spin-dependent hoppings. 

II. BAND-STRUCTURE ENGINEERING 

To set the stage, let us first cast the analysis of 
graphene's honeycomb lattice in a form that lends it- 
self to generalization. The bipartite honeycomb lattice 
has two atoms (A and B) per unit cell; each sublattice 
forms a triangular lattice and the hopping Hamiltonian 
in Fourier space takes the form 

H-\ */(k)l .2) 

where t is the hopping amplitude and a the intercarbon 
distance, while /(k) = 1 -I- 2 exp(i3/cya/2) cos(-\/3fca;a/2). 

At half-filling the Fermi surface consists of two 
inequivalent Dirac points K and K' at momenta 



±(27r/3-\/3a, 27r/3a) = k±, respectively. Expanding 
around these points leads to two copies of the two- 
dimensional Dirac equation for S = 1/2 with the sub- 
lattice providing the (pseudo)spin degree of freedom: 

^""2(10)' ^^=2(^0) 

vp = 3ta, the missing factor 1/2 as opposed to graphene^ 
arising since the eigenvalues of the spin are ±1/2. 

If we now think of the honeycomb lattice as a layered 
structure, with the A and B triangular sublattices off- 
set in height by an amount h, it is natural to ask what 
happens if one adds a third, and then further, triangular 
layers (see Fig. [T|). For small h, the nearest neighbours 
of a given site are in the layers directly above and below. 
When the layers are stacked in the sequence of the face- 
centred cubic lattice in a [111] direction, one obtains a 
band diagonal hopping Hamiltonian for a system of 25* -1-1 
layers: 



Hs = t 







aoi/(k) 



a*J*ik) 
alJ*{k)^ 






ai2/(k)^ 

0^ a2S-i,2s/(k) 
45-i,2s/*(k) 



(4) 



where we have allowed for different interlayer hopping 
strengths by introducing the q's. Indeed, regardless of 
the choice of Ofi^i+i, several properties of the spectrum of 
Hs immediately follow from the form of its characteristic 
polynomial Cs{X, k) = det {Hs — Al) which reads: 



Cs(A,k) = \tf {k)f^+^ det (a+ + Q-) - Al 



(5) 



Firstly, near the Dirac 



points {K and K'), all bands are linearly dispersing, 
simply because |/(k)| oc |k — k±|. Secondly, for inte- 
ger S the Hamiltonian must display a flat band. The 
matrix a"*" + a~ possesses a symmetric spectrum (if A 
is an eigenvalue, then so is also —A). The Hamiltonian 
also possesses this symmetry, which can be phrased as a 
chiral symmetry I]iJ(k)I]t — — i/(k) with S a unitary 
operator, as further discussed in the appendix [A] For 
an odd number 25" -|- 1 of bands the chiral symmetry im- 
plies that (at least) one eigenvalue A must be zero, which 
translates into a flat band A = |/(fc)|A = 0.— 

If now, in addition, we choose the interplane hopping 
amplitudes so that a+ ~ S^, where = -\- iSy is 
the raising operator for spin 5, we obtain a spectrum 
Enik) = nt\f{k)\, where n = -S,-S + 1,...,S. This 
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requires placing the adjacent layers at certain distances 
from each other, so that the overlap of the wavefunctions 
would produce the appropriate hopping integrals between 
subsequent layers, whose relative strength is further spec- 
ified in Eq. ©. 

As a result, we obtain an effective Hamiltonian near 
the K point. 



^^s(p) = vpS ■ p. 



(6) 



where A = A/|/(k)| is independent of k, and so is Here 



p = k — k_|_ measures the (small) distance from the 



Dirac point at K, and similarly for the K' point. If we 
finally add a potential of strength ASz, which can in prin- 
ciple be generated straightforwardly via an electric field 
applied perpendicular to the layers, representing distinct 
chemical potentials for each layer, we have 



Enip) 



(7) 



For integer S*, n = invariably corresponds to a flat band, 
however, no longer due to chiral symmetry but due to a 
less general symmetry that is specific to the low-energy 
Dirac-like Hamiltonians and that requires fine-tuning of 
the parameters Uij. Fortunately for experimental real- 
izations, rather natural setups like equidistant layers will 
satisfy the conditions for a fiat band, as discussed in the 
appendix 1X1 in more detail. 
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FIG. 1. (Color online) The schematic representation of the 
family of lattice models leading to the spin-S Dirac equation 
is shown from the side (a) and from above (b). The dashed 
lines denote interlayer hopping processes, while the intralayer 
thin solid lines are guide to the eye, emphasizing the pla- 
nar triangular structure, but do not represent any hoppings. 
The lowest and highest (blue A) planes are exactly on top of 
each other. For S = 1/2, only two adjacent layers need to 
be considered (e.g. red A, and green B), for S — 1, three 
neighbouring layers (e.g. lower blue A, red B and green C), 
for S = 3/2, all four layers, while for higher 5's, one needs to 
continue up- or downwards with ABCABC... stackings. Note 
that the interlayer hoppings should be unequal to realise the 
perfect Dirac equation, Eq. ^ 



It can now be verified straightforwardly that the cases 
5 = 1/2 and S = 1 correspond to the known instances of 
the honeycomb and dice (or T3) lattices^i^, respectively. 
The S* = 3/2 case for four layers reads 

V3i/(k) 

x/3tr(k) 2tf{k) 

2t/*(k) VStfik) ' ^ ' 

\/3i/*(k) 

which can be supplemented with an additional gap, com- 
ing from AS:, with = diag(3/2, 1/2, -1/2, -3/2). A 
very similar lattice structure has been proposed in Ref. 

Note that the simple form of the Hamiltonian also gives 
immediate access to the wavcfunctions in layer (pseu- 



dospin) space, as its eigenfunctions are obtained from a 
simple rotation in spin space: the quantization axis of S 
is given by an effective field direction, h, whose compo- 
nents are given by Re/(k), Im/(k) and A, respectively. 
The various a-prefactors above and below the diagonal 
arc chosen according to the conventional matrix repre- 
sentation of the spin matriceai^. For example, above the 
diagonal, the matrix elements of the raising ladder oper- 
ator appear as 

{n'\S+\n) = Sn'.n+iVS{S + l)-n{n + l), (9) 

where Sz\n) — n\n). The resulting spectrum consists 
of equidistant energy levels at each given momentum. 
As mentioned above, with a different choice of a, the 
spectrum would still be linear. The Dirac cones are ro- 
bust in this sense. Furthermore, although the bands at 
a given momentum would not necessarily be equidistant 
anymore, the flat band will still survive if certain sym- 
metries are present, as discussed in the appendix \X\ 

For A = the wavefunction corresponding to the flat 
band is such that the probability of finding a particle in 
even layers is exactly zero. For example, in an 5 = 1- 
trilayer, the red plane, sandwiched between the blue and 
green ones (see Fig. [T]), is completely blocked for the flat 
band wavefunction. 

The stacking we propose here is of course quite famil- 
iar. A succession of triangular planes ABCABC... as 
displayed in Fig. [U is just the face-centred cubic lat- 
tice viewed along a [111] direction. Another stacking, 
ABABA..., corresponds to the hexagonal close packed 
lattice structure. The hopping Hamiltonian, Eq. is 
simply modified to take into account this stacking: /(k) 
is replaced by its complex conjugate for hopping BA, CB 
or AC. For example the stacking ABCB would result in 

aoi/(k) - 

^ a*J*ik) ai2/(k) 

' alJ*{k) a^3r(k) ' 

a23/(k) . 

(10) 

We will say that the chirality between the third and the 
fourth layer has been flipped. 

The spectrum (and in particular Dirac cones and flat 
bands) is not affected by this change, which affects only 
the phase of the matrix elements. Indeed, a unitary 
transformation, changing the i*^ spinor entry ipi{k) — > 
exp[2i arg/(fc)]V'i(fc) changes f*{k) f{k) in 

Around half-filling, where the continuum description 
applies, this corresponds to fiipping the chirality in the 
Hamiltonian between adjacent layers. This change of chi- 
rality is at the origin of a change in the Berry curvature 
(detailed in Sec. |V|. Note that one can successively 'fix' 
the phases of the off-diagonal terms to agree with a refer- 
ence stacking without altering diagonal terms which may 
be present. 
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Finally, we emphasize again that for a layer separation 
h < a-\/2/3, where a is the triangular lattice constant, 
our Eq. Q represents nearest-neighbour hoppings only. 

Having outlined a path towards general lattices with 
Dirac physics, we next discuss some of the basic proper- 
ties of such electronic systems. 



III. DENSITY OF STATES 

The density of states (DOS) for A = is given, using 
the low energy Dirac Hamiltonians, by 



p,n— — 5 



n>0 



(11) 



per spin, valley and unit cell, Ac being the unit cell area. 
The DOS remains linear in energy, similar to graphene, 
but exhibits a sharp peak due to the flat band^ii^ for 
integer spin realizations. The DOS can be simplified to 



2tt vi 



6 



Sic,) 



(12) 



for integer spins, and 



Ac \U1\ I TT' 



(13) 



for half- integer spins, where "^{x) is Euler's digamma 
function. Due to the momentum integral in Eq. (|lip . 
these results are only valid for |fc| <C kc with kc the cut- 
off, which translates to ^ vpkc- Note that, for large 
spin 5* ^ 1, the maximal slope of the DOS right at the 
Fermi energy is 7r^/8 times larger than for spin- 1/2 for 
half integer spins, and 7r^/6 times larger than for spin- 
1 for integer spins. The S dependence of the slope of 
the DOS is shown in Fig. which changes very little 
with S in the integer or half-integer sector. With in- 
creasing 5, additional Dirac cones appear with increas- 
ing slope, thus with a much reduced contribution to the 
DOS. As opposed to that, these high energy bands con- 
tribute more at high energies, since their bandwidth also 
increases with S. 

The original lattice model provides us with additional 
features, not captured by the low energy approximation, 
such as the presence of van Hove singularities around 
uj = nt with n = —5, —5 + 1 . . .S but n 7^ as p{ui) ~ 
ln(aj/|n|t). In addition to the two peaks for graphene 
with S ~ 1/2, increasing number of additional pairs of 
peaks appear in the DOS for 5 > 1. 
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FIG. 2. (Color online) The slope of the linear in energy den- 
sity of states is plotted for various values of S, the dashed 
black lines denote the asymptotic values, 7r^/6 and 7r^/2 for 
integer/half-integer spins, respectively. 



IV. OPTICAL CONDUCTIVITY 

Another characteristic quantity of Dirac fermions is 
the optical conductivity, which, for S = 1/2 at half fill- 
ing and r = 0, is completely structureless and constant. 
In the presence of additional Dirac bands, new interband 
transitions occur. The current operator in the x direc- 
tions is given by jx = vpSx- Its equation of motion as 
well as those of the other spin components are 



dtSx = VFPySz, 

dtSy = -vpPxSz, 
dtSz = vp [pxSy -pySx) 



(14) 
(15) 
(16) 



for a given momentum. This is easily solved for Sx (t) as 



Sx{t) = Sx [sm^{Lpp) cos{v Fpt) -f cos^((^p)] -|- 

+ -5*^ sin(2</?p) [1 — cos{vppt)] + sin((y9p) sm{vppt), 

(17) 

where tan((^p) = px/py 

The current-current correlation function is evaluated 
from this as 



Tl,p 

= 2i E]('S'j,) sm{ipp) sm{vFpt) = 

n,p 

— 2i''^^npsm'^(ipp) ai'a{vppt). (18) 



After Fourier transformation, the optical conductivity 
contains two parts as 



(19) 



5 



per electron spin and valley, and the Drudc weight is 



D^l^L5' + l/2jlnf2cosh^ /i 



2kBT 



(20) 



which agrees with that of graphenei^ for S = 1/2, while 
the interband part reads as 



9 

en 



s 

n=-S 



{21] 



where f{x) = l/(cxp((.T — fi)/kBT) + 1) is the Fermi 
function, n the chemical potential, [a;J denotes the in- 
teger part. Since the particles residing on the flat band 
cannot propagate, their group velocity is zero, so that 
their contribution vanishes to the Drude weight. This 
explains the integer part function. On the other hand, 
they have a finite matrix element between adjacent levels, 
and contribute to interband transport, which contains all 
allowed 25* processes between 25" + 1 levels. 

At the Dirac point (/i = 0) at T = 0, the Drude weight 
disappears, and the interband conductivity reads 



Winter (^) 



l/T I 2 



- half-integer S 
8 ^ 
integer S 



(22) 



Away from the Dirac point, \_S + 1/2J interband transi- 
tions are allowed, as can be checked in Figs. [3]and|4l 
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FIG. 3. (Color online) The interband part of the opti- 
cal conductivity for the spin-5 Dirac equation is shown for 
ksT / ^ = 0.0125 for several values of S. The number of pos- 
sible interband transition is [S + 1/2J. 

The calculated intra- and interband optical conductiv- 
ities differ significantly from those in graphene. First, 
the interband part is sensitive to the number of bands 
and away from half-filling, several steps are possible as 
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FIG. 4. (Color online) The band structure of 5 = 3/2, visual- 
izing the minimal frequency of the allowed optical transition. 
As opposed to graphene, transitions well below u = 2^ are 
possible. 



opposed to graphene, where only a single step is allowed. 
Second, the universal value for the optical conductivity 
at half filling and finite frequencies is proportional to 
pseudospin-5 value, which should also affect the trans- 
parency as 



r= 1 



1 — naQED 



5(5 + 1) + 



27: 

— cr(w) 
c 

( 1 



half-integer 5 
integer 5 



, (23) 



where the lower line is obtained upon Taylor expanding 
the upper line, and is only valid for 5 < 2. Here, val- 
ley and physical spin degeneracies are included, uqed — 
/he is the fine structure constant, c the speed of light. 
For 5 = 1/2, this reproduces the T — 97.7% optical 
transparency of graphene. Therefore, the universal value 
of the optical response immediately reveals the underly- 
ing pseudospin-5 structure, as shown in Fig. [SJ Third, 
while the interband response takes the contribution of 
the flat band into account, the intraband one (Drude) is 
insensitive to its presence due to the zero group velocity 
of the fiat band. 



V. TOPOLOGICAL PROPERTIES OF THE 
BAND STRUCTURE 

Here wc discuss the topological properties our model 
in the absence of a gauge field. It will turn out that 
there can be topologically non-trivial grounds states. 
The topological invariant we study — the spin Chern 
number — depends on the number of layers and on the 
band fillings. However, it does not appear to depend on 
the stacking configuration, although the Berry curvature, 
from which it is calculated, docs. 
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FIG. 5. (Color online) The universal optical transparency is 
shown for half filling as a function of S. The blue circles 
denote the exact expression from Eq. (|23p (upper line), while 
the black squares come the approximate formula, valid for 
smaU 5" (lower line in Eq. fl^ ). 



We study the integral of the Berry curvature, where the 
contribution to the Berry curvature from a given band is 
given by 

^ ^ J2 [{n, k\dk^H\n', k) (n\ k\d,^H\n, k)] 

n=-S, 
n 



[{En{k)-En'{k)Y 



(24) 



n is a band index and k) a single particle eigenstate of 
H with eigenvalue £'„(k). In the thermodynamic limit 
the summation over momentum turns into an integral. 
When this integral goes over a compact manifold like the 
Brillouin zone, one obtains a topological invariant called 
the first Chern number 



BZ 



d^. 



-B„ 



(25) 



A non-zero Chern number requires breaking of time 
reversal symmetry (TRS). This can be accomplished by 
an external magnetic field like in the integer quantum 
Hall effect^, but also by gap terms that break TRSi^. 
Settings of the latter kind with flat bands with non-zero 
Chern number have been reportedi^r— . In this section 
we only study gap terms that preserve TRS, hence the 
Chern number is zero. However, the Chern number is in- 
tegrated from a non-trivial Berry curvature that derives 
mainly from the two singularities in the Brillouin zone - 
the Dirac cones at the K and K' points. When calcu- 
lating the total vorticity of a configuration of quantized 
vortices, it is usually enough to calculate the vorticity of 
the individual vortices as if they were isolated, and then 
add up the quanta (including the sign) to obtain the to- 
tal vorticity. Similarly, it is usually enough to calculate 
the integral of the Berry curvature for individual Dirac 
cones, living not on a Brillouin zone but on the infinite 
plane of momenta and then add up their contributions 
to obtain the same result which would originate from in- 
tegrating the full band structure over the Brillouin zone. 



We will use both approaches and demonstrate explicitly 
(in Table H]) that they match. 

Another topological invariant is the spin Chern num- 
ber. It is calculated from the contributions to the Chern 
number individually for the two components of the phys- 
ical spin, but instead of adding up the two contributions 
to get the Chern number, one takes the difference to ob- 
tain the spin Chern number. 

The charge Chern number is related to a topologically 
quantized Hall current of charge^^. 



(26) 



with the summation taken over filled bands. Likewise, 
the spin Chern number can be used to determine the 
quantized spin-Hall conductivity. This is the case with 
the intrinsic spin-orbit coupling (SOC) of graphcne^^ 



HsO = ^TzSzUz 



(27) 



which preserves TRS. Here and cr^ refers to the valley 
and physical spin degrees of freedom. We will generalize 
this SOC to arbitrary pseudospin-S*. 

In addition, the different stacking configurations lead 
to interesting changes in the integral of the Berry cur- 
vature in the single cone approximation, which calls for 
future research on the role of stacking order on topolog- 
ical properties. 



A. Single pseudospin-S Dirac cone 

Let us first focus on the integral of the Berry curvature 
of a single spin-S* Dirac equation for spinless electrons. 
In itself, it is usually not a topological invariant, but the 
topological invariants can often be understood in terms 
of the contributions for the single cones, and the latter 
can in the simplest case be evaluated analytically. In 
some cases, these non-zero contributions add up to zero, 
as must be the case for the Chern number in the time 
reversal symmetric setup that we consider. However, we 
will also see that for the spin Chern number, another 
topological invariant, the contributions add up to an in- 
variant that can be nonzero. 

For a given band in the spin-5 Dirac equation, assum- 
ing a gap of the form AS'^, the integral of the Berry 
curvature is evaluated around the K point by assuming 
an isolated Dirac cone in the momentum plane in the 
expression of Eq. (f24| . Thus 



K 



nAvppdp 



[{vFpy 



l3/2 



nsign(A), (28) 



where n takes the allowed values of Sz, which also in- 



dexes the bands. For the S 



Dirac equation, this 



reproduces the known result^ C± = ±isign(A) with 
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the upper/lower sign corresponding to the upper/lower 
Dirac cone. 

The contribution to the transverse conductivity from 
an individual cone at the Dirac point K is the sum of 
the above C^'s from the filled bands. In the Dirac cone 
approximation a gap exists only around zero energy (be- 
tween n = —1/2 and n = 1/2 for half-integer spin or 
between the bands n ~ —1, n = and n = 1 for inte- 
ger spin). Therefore, as long as the chemical potential 
satisfies \fi\ < A/2 for half- integer and < \fj,\ < |A| for 
integer spins, one obtains also in the low-energy approxi- 
mation a half-integer quantized transverse response (per 
spin and valley) 



ri<0 



K 



' h 2 
for integer S 



sign(A)- 



- ^ < signj^) half-integer S 



(29) 



An increasing chemical potential will cut into some bands 
which destroys the half- integer quantization of a^y . 

By choosing a stacking different from ABCAB..., other 
decompositions across the bands can be obtained for the 
Berry curvature in the single cone approximation, which 
lead to a result different from Eq. (|28p . For example, 
for 5' = 1, i.e. the trilayer with ABC stacking, we obtain 
for positive A that (Cf ,C^,C5i) = (1,0,-1) per spin 
and valley. By flipping the chirality between the second 
and third layer (which corresponds to an ABA stacking), 
we obtain instead (1/2, —1, 1/2) as integrals of the Berry 
curvature for the successive bands. (These integrals were 
evaluated numerically, in contrast to those related to the 
spin-S" Hamiltonians.) Most importantly, the contribu- 
tion of the fiat band gets modified to 1. Various stacking 
patterns, whose variety grows with the number of layers, 
are listed in Table H] As we demonstrate next, the inte- 
gral of the Berry curvature around a single cone can be 
used to determine the spin Chern number of the lattice 
model. 



B. Chern numbers on the lattice 

To obtain a topological invariant the full Brillouin zone 
(including the K' point) has to be considered. For a set- 
ting with TRS, which we consider here, the Chern num- 
ber has to be zero, even though the Berry curvature con- 
tributions around the individual cones may be nonzero, in 
which case they cancel. This cancellation can be avoided 
by involving the physical spin: by generalizing the intrin- 
sic SOC of graphcne in Eq. (|27l) to general spin-S". we 
can end up with non-zero spin Chern number and thus 
finite spin Hall conductivity. 

By going back to the original lattice model and defining 
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(1/2,-1/2) 


(1,-1) 


1 


ABC 


(1,0,-1) 


(2,0,-2) 


1 


ABA 


(1/2,-1,1/2) 


(0,0,0) 


3/2 


ABCA 


(3/2,1/2,-1/2,-3/2) 


(3,1,-1,-3) 


3/2 


ABAB 


(1/2,-1/2,1/2,-1/2) 


(1,-1,1,-1) 


3/2 


ABCB 


(5/4,-1/4,-5/4,1/4) 


(1,1,-1,-1) 


2 


ABCAB 


(2,1,0,-1,-2) 


(4,2,0,-2,-4) 


2 


ABCAC 


(15/8,1/2,-3/4,-3/2,-1/8) 


(2,2,0,-2,-2) 


2 


ABCBC 


(11/8,0,-3/4,0,-5/8) 


(2,0,0,0,-2) 


2 


ABCBA 


(5/4,-1/2,-3/2,-1/2,5/4) 


(0,0,0,0,0) 


2 


ABABA 


(1/2,-1/2,0,-1/2,1/2) 


(0,0,0,0,0) 


2 


ABACA 


(0,-1,0,1,0) 


(0,-2,0,2,0) 



TABLE I. The integral of the Berry curvature for a single 
pseudospin-5' Dirac cone (C^) and the spin dependent Chern 
number (Cj^) for the lattice model are shown for all possible, 
non-equivalent stacking patterns for S < 5/2 and A > 0. 
These are evaluated numerically for both the continuum and 
lattice model, and their relation can be checked using Eq. 
(I33ll. Note that C2 = -Ci. 



the full lattice version of Eq. ([27)1 following Ref. [22 as 
2A 



so .lattice 



3V3 



2 sin{\/3kxa/2) coa{3kya/2) 



(30) 



which should be added to Eq. (j4]) , we get the spin Chern 
numbers, from Eq. (|24|) (with a numerical integration 
over the entire Brillouin zone) as 



C,t = C^-Ci-4nsign(A) 



(31) 



where the factor 4 comes from the valley and physical 
spin degrees of freedom, and C is the Chern number for 
up (cr =t) or down {a —l) spins and CJ^ = — C^. 

Whether the system qualifies as a spin-Hall insula- 
tor is decided^^ by the Z2 invariant ly, defined by ly = 
^„ C^/2 (mod2) = ^„ 2nsignA (mod2), and summa- 
tion is taken over filled bands. 

We thus find that integer pseudo-spins contribute with 
even numbers to the sum in i/, and are topologically triv- 
ial for the lattices in Fig. [U giving = 0. By con- 
trast, half-integers spins contribute with odd numbers to 
the sum, and having an even number of filled bands at 
half filling (i.e. S = 3/2, 7/2, 11/2. . . ) adds up to an 
even number, thus again v ^ 0. As opposed to this, 
half-integer spins with 5=1/2, 5/2, 9/2... have an odd 
number of filled bands at half filling, adding up to an odd 
number, resulting in v = 1 and topologically non-trivial 
behaviour. 

The spin dependent Chern numbers of the lattice 
model can be obtained from the single cone results of 
the previous section. In the single cone approximation 
the spin-orbit coupling of Eq. ([50]) simplifies to Eq. (|27l) . 
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These single cone contributions to the Chern number sat- 
isfy 

(A) = -Cf (A) = C^„(-A) = -Cii-A). (32) 

The opposite sign of the gap term for the two Dirac points 
derives from in Eq. (j27|) . 

By taking both Dirac points into account, we obtain 
from the single cone results the spin Chern number 

d = C^{A) + Cf (-A) = C,f (A) - C^^JA). (33) 

The result agrees with the result in Eq. (pij) found from 
integrating over the Brillouin zone using the full band 
structure. Eq. (|33p immediately implies that the spin 
Chern number of the flat band is zero (Cq = 0), regard- 
less of the value of C^{A). 

We can also consider other stacking patterns, as we did 
in the single cone case. The correspondences in Eq. ((33|) 
hold for arbitrary stacking patterns on the lattice. For 
example, the trilayer with ABA stacking with A > 
yields zero spin Chern numbers for all bands, unlike the 
S = 1 case derived from the ABC stacking. However, like 
the ABC stacking, the ABA gives a spin Chern number 
that is topologically trivial. Results for other stackings 
are shown in Table HI One interesting general conclusion 
that we can draw based on these results is that while the 
single cone Chern number contributions for each band 
are redistributed significantly with different stacking pat- 
terns, this does not affect the Z2 topological invariant. 
This invariant will therefore be determined only by the 
number of layers but not by the stacking. We have also 
checked that all non-equivalent stackings for the S* = 5/2 
case (not shown here) give 1/ = 1 at half filling. The 
invariance to changes in stacking applies also away from 
half-filling, as long as the chemical potential lies between 
the bands. 

We can also engineer nearly flat bands with non- 
trivial topology, similarly to Refs. [T7l - [l9t when A S> 
all bands become practically fiat as En{p) ~ nA + 
nt^|/(k)p/2A, i.e. the hopping occurs only to second 
order in perturbation theory. Therefore, it becomes pos- 
sible to fill the separate bands one by one. Then, for ex- 
ample, the quarter filled S ^ 3/2 case, which corresponds 
to a completely filled E_^/2{p) band, becomes topologi- 
cally non-trivial with i/ = 1 for all stackings. Note that 
when A <C t, quarter filling in this case gives partially 
filled £'-3/2(p) and E_i/2{p) bands. In the same vein, the 
1/3 filled 5 5/2 lattice with flattened bands (A > t) 
is topologically trivial with v — 0. Thus, a trivial ground 
state can become non-trivial (and vice versa) when the 
chemical potential is lowered or increased to the next 
band gap. 

Another observation we have made is that while the 
topological invariants (the Chern number and the spin 
Chern number) are robust with respect to the variations 
of a's in Eq. the integrals of the Berry curvature for 
a single cone (C^) are not invariant. However, for some 
stacking patterns, is rather insensitive to changes in 



a's. In particular, the single cone results for ABABA 
and ABACA stackings are also recovered for uniform in- 
terlayer hoppings. 

We close this section with the remark that topological 
invariants do not depend only on S but also on the num- 
ber of non-equivalent Dirac cones and thus on the specific 
form of the lattice. For example, T3 and Lieb lattice a^'^i^^ 
with even and odd number of 5* = 1 cones, respectively, 
belong to different Z2 class^. In the presence of intrinsic 
SOC, the T3 lattice with two inequivalent cones, pos- 
sesses the trivial Z2 index. On the contrary, the Lieb 
lattice has a single cone in its band structure and has 
therefore a ground state with a non-trivial Z2 invariant, 
and realizes a spin-Hall insulator in the presence of SOC. 



VI. TOPOLOGICAL PROPERTIES IN THE 
PRESENCE OF A MAGNETIC FIELD 

A topological property of our lattice that turns out to 
depend dramatically on the stacking configuration is the 
number of zero-modes in a magnetic field. In the case of 
a uniform magnetic field, these zero-modes are nothing 
but E = Landau level (LL) states. We now show that 
by changing the stacking from the ABCAB... stacking to 
some other stacking, the E = Landau level degeneracy 
will increase by a multiple. 

This can immediately be seen in experiment as an in- 
creased step at = in the steps of quantized Hall con- 
ductivity axy as a function of chemical potential, in a 
way analogous to the well-known examples of mono- and 
bilayer graphene: in the former, all Landau levels have 
the same degeneracy, while in the latter, the degeneracy 
of the £' = LL is twice the one of the other o^^i^^ . Such 
degeneracies and their lifting play an important role, for 
example at integer fillings in the context of multicompo- 
nent quantum Hall ferromagnetism, see e.g. Ref. [27l 

The reason for this increased multiplicity thanks to 
restacking is rather easy to understand. Consider the 
gauge invariant momentum operator I?i(x) = —idx^ — 
Ai{'x). The existence of zero-modes relies on the fact one 
of the chiral Dirac operators 2?± — ~i{dx ± idy) — {A^ ± 
iAy) has a non-trivial kernel when the net flux of A is 
bigger than one flux quantum. The number of states in 
the kernel is given by the number of flux quanta. De- 
pending on the sign of the total flux one finds non-trivial 
solutions either to V^ij} = or to 2?_-i/' — 0, but not to 
both. This latter fact comes into play in an interesting 
way when we start to flip the chiralities of the hoppings 
between individual layers by restacking, as we will now 
come to. 

Like in graphene we discuss the Landau level spectrum 
in the linearized regime of the low-energy Hamiltonians, 
that is, in terms of the Dirac cones. (In the case of the 
full band structure we cannot write the Hamiltonian only 
in terms of the chiral combinations 'D±, which is neces- 
sary for the analytical discussion of zero-modes.) The 
Dirac Hamiltonian is then written in real space and the 
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magnetic field is introduced by minimal coupling. For 
the ABCA stacking, we find at the Dirac point K 



Ha 



BCA 



















"23^+ 










(34) 



With the ABAB stacking we instead have 



Ha 



BAB 







"01^4 



aoiV^ 
ai2'D+ 



'12 


















, (35) 



that is, with the chiralitics of the matrix elements relating 
the second and third layer flipped. Thus, by restacking it 
is possible to obtain several columns with only one chi- 
rality of Dirac operators and not both. Such columns 
will contribute with new zero-mode solutions and will in- 
crease the zero-mode degeneracy by a multiplicity factor. 
Assume that the flux is such that there are n solutions 
ipi {i = f,...,7T,) to V^ip = and hence no solutions 
to V-ip = 0. Then there are only the n zero- modes for 
^^ABCA in Eq. ([M]) of the form 



* = (V., 0,0,0)^ 



(36) 



^^ABAB in Eq. p4[) . on the other hand, has the same n 
zero-modes, but also n additional zero-modes of the form 



vI/ = (0,0,V.,0)^ 



(37) 



which are not solutions to -//abca because of the mixed 
occurrence of T)^ and I?_ in the third column of i?ABCA ■ 
Thus, the Hamiltonian i?ABCA has twice as many zero 
modes. For larger number of layers one has an even big- 
ger number of different stacking configurations to choose 
between, each with different implications for the zero- 
mode degeneracy. The one extreme case is given by the 
"face-centred cubic" stacking ABCABC... (a configura- 
tion without flipped chiralitics), which remains at the n 
zero-modes for an arbitrary number of layers. The other 
extreme is given by the hexagonal close-packed stack- 
ing ABABAB... (with an alternating sequence of chirali- 
tics), where the number of zero modes of 25 -I- 1 bands is 
n{S+l/2) for half-integer S, and nS or n{S+l) (depend- 
ing on the sign of the flux) for integer S. Other stack- 
ings give some intermediate multiple of n zero-modes. 
A quick inspection shows that the contribution from the 
other Dirac point just duplicates this result for any stack- 
ing configuration, thus there will be a factor two due to 
valley degeneracy. 

The above observed flexibility to increase the zero- 
mode degeneracy by a simple change of stacking is in 
stark contrast with the situation in multilayer graphene. 
Multilayer graphene has in the simplest approximations 
indeed multiple times the degeneracy of the monolayer 
Hamiltonian^^. However, although the structure of the 



Hamiltonian depends sensitively on the stacking^S, the 
degeneracy turns out in the simplest approximation to 
be independent of stacking, even in the case of a non- 
uniform vector potentiali^ 

These results in fact also apply for magnetic fields 
which are no longer uniform. While in this case. Landau 
level degeneracies will in general be lifted, the E = Lan- 
dau level for 2d Dirac electrons is an exception. The rea- 
son is the widely known general property of Dirac opera- 
tors in a vector potential of arbitrary distribution which 
is treated by the Atiyah-Singer index theorero^, but can 
be understood also on a less formal level thanks to the 
neat argument by Aharonov and Casher— . This may 
imply a £' = Landau level that is qualitatively sharper 
than the other Landau levels, since only the latter are 
broadened by the non-uniform component of the mag- 
netic field. In graphene a non-uniform component can 
be due to the effective magnetic field introduced by the 
corrugation of the graphene membrane. That the E = 
level remains relatively sharp in graphene has been ob- 
served in experimentsi^'^. Even already for weak mag- 
netic fields, one can expect an increased density of states 
at i? = due to these zero-modes. If the degeneracy of 
the zero-modes can be multiplied, as we have shown for 
our example, then such a peak in the density of states 
should grow with the same multiplicity. 

Notice that the Landau level degeneracy depends 
strongly on the chosen stacking of the layers, thus also 
influencing the height of the zero energy peak in the den- 
sity of states. This is in stark contrast to the zero field 
results, where the DOS is stacking independent, as stud- 
ied in Sec. HI. 



VII. KLEIN TUNNELING ON A POTENTIAL 
STEP 

Finally, we discuss Klein-tunneling of spin-S* Dirac 
electrons. As we have seen, qualitative differences arise 
between half-integer and integer spins. The transmission 
amplitude for spin- 1/2 Dirac electrons has been studied, 
in connection to graphene, in Ref. [s^. The spin-1 case 
and the influence of the flat band was studied in Refs. 
I^,!!!, and[36l and all-angle perfect transmission was found 
at specific energies. Here we discuss Klein tunneling for 
the pseudospin-1 case, and focus on tunneling into the 
flat band, as shown in Fig. [6] 

But before doing so, let us discuss the general matching 
conditions of the wavefunction for general pseudospin- 
S. In this case, the spinor wavefunction takes the form 
(omitting spatial coordinates for simplicity) 



^ = (*l,^'2,...,^'2S,*2S+l) 



(38) 



To determine the matching conditions^, we integrate the 
eigenvalue equation of the Hamiltonian in Eq. or 
more generally from Eq. (U), after expanding it around 
the K point, Hs'i! — E'i/ from x = —xq to x — xq and 
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send xq to zero. As a result, we obtain, assuming non- 
diverging scalar and vector potentials 

*2(-:co) -*2(a;o), (39) 
a*_i_j*j(-a:o) + aj,i+i'I'j+2(-2:o) = 
= a*_i ,^z{xo) + Q;i,j+i*i+2(a;o) for i = 1 ... 25 - 1, 

(40) 

*2s(-2:o) = *2s(a;o)- (41) 

In general, this implies that \I/2 and \I'2S must be contin- 
uous, since there is only a single non-zero element in the 
first and last row of the Hamiltonian matrix. In addition, 
Eq. (|40p . involving other components of the spinor, con- 
tains only and \I'i+2 , hence there is no mixing between 
even and odd components. 

The case of half-integer pseudospin (with even 2S +1) 
implies that the two continuous spinor components, Eqs. 
([M)) and (|¥T|) contain one even (2) and one odd {2S) in- 
dex. Therefore, the continuity of e.g. 4*2 implies through 
Eq. (|40|) that of ^'4. The continuity of 4*4 in turn implies 
the continuity of ^'e and so on. The very same procedure 
can be carried out for the odd components. Therefore, 
each component of the spinor changes continuously, thus 
the whole wavefunction remains continuous across a po- 
tential barrier. 

The case of integer pseudospins is different: 25 -f 1 is 
odd, therefore only two even components (\I'2 and ^'25) 
arc required to change continuously (as opposed to the 
one even and one odd components for half-integer S). 
This implies that all even components must be contin- 
uous, but there are only S equations for the remaining 
5-1-1 odd components. Only the continuity of the linear 
combinations of neighbouring odd components, Eq. (PO)) 
with i = 1, 3,. ..25 — 1 is required across a barrier, but 
nothing can be said about the individual components. 

Similar considerations apply along the y directions, in 
which case the wavefunction still changes continuously 
for half-integer 5, while only the even components remain 
continuous for integer 5. Implicitly, this difference can 
be traced back to the absence or presence of a flat-band. 

Here we consider the scattering of pseudospin- 1 elec- 
trons on a sharp potential step of the form VQ{x), V > 0. 
The case when the energy of the injected electron dif- 
fers from V has already been considered in Refs. H and 
IsgI . However, when the energy of the incident electron 
is exactly E = V, scattering into the non-propagating 
flat-band becomes possible. The electrons on the flat 
band do not possess a well defined Fermi surface, since 
all particles residing on the flat band have identically zero 
energy. Thus, an incident electron with E = V can be 
scattered to any momentum state of the flat band within 
the barrier. 

At normal incidence {ky = 0), there is perfect trans- 
mission (T = 1), since transmission to the upper or lower 
Dirac cones is possible (excluding the flat band). 

At a finite angle, scattering to the propagating cones 
is forbidden by momentum conservation (ky does not 
change). In this special case, the wavefunction on the 




FIG. 6. (Color onlme) A sharp potential barrier for the spin-1 
Dirac equation, the thick red lines denote the non-propagating 
zero energy states, while the short black arrows stand for the 
velocity of the branches. 



left and right hand side of the barrier (suppressing the 
cxp{ikyy) term) is given by 



*^(.T < 0) = - 



^^{x > 0) = y ^ 



' exp(j(/3A;) 

V2a I exp(ifca;x)-|- 
\ exp(-i(pfc) 

-cxp(-iv3fe) \ 

\f2oL cxp{-ik^x), (42) 

-exp(i.^fc) J 

exp(ifc^a;)-|- 
-exp(-iy3fc/) / 



e{ky) 

a \ I cxp(-|A:j^|a;) 

e(-fc,) 



(43) 



where ky is conserved, i.e. ky 



k'y, V 



v\k\, tan((pfc) 



kx/ky, tan((y5fc') = k'^/ky, and the lack of Fermi surface 
implies that any state on the flat band is available for 
transmission without any restrictions on fc^, explaining 
the summation over k'^, a ~ 1. The last term describes 
an evanescent mode in the flat band. Applying the con- 
tinuity of ^'2 from Eq. ([59]) implies that the reflection 
coefficient r = — 1, from which the reflection probability 
is R= \r^\ = 1, and r = 0. 

As far as such a stationary solution is concerned, states 
in the flat band to the right of the barrier may also be oc- 
cupied. However, as the group velocity on the flat-band 
is zero, the transmission probability is also zero, or in 
other words, the probability current is zero through the 
barrier. The resulting picture thus consists of standing, 
non-propagating waves, extending to both sides of the 
barrier: on the left, it is made of two counterpropagat- 
ing waves (in the x direction), whose interference leads 
to a standing wave, while on the right, the zero energy 
mode is non-propagating by its very nature. Although 
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the wavefunction to the left of the barrier is uniquely 
determined by specifying energy E = V and perpendic- 
ular momentum ky, its corresponding part to the right 
of the barrier has many degenerate versions due to the 
flat band, and can thus host a large number of different 
states ^ L, even for a fixed ky. 

This can be made explicit as follows. Using Eq. (|40)) 
to connect Eqs. (|42|) and (j43| . we get 

cos(^fc) = ^ + -J=^t(fc')sin(v5fc'), (44) 

which can have ^ L distinct set of independent solutions 
in terms of {t{k')} and a. Crucially, each such solutions 
corresponds to zero transmission probability and perfect 
reflection. 

Very similar considerations apply to the case of -E = 0, 
namely an electron in the flat-band to the left of the 
barrier, scattered to propagating states to the right. 

In the presence of many bands (S* > 1), tunneling be- 
tween them occurs with a greater variety, and interband 
tunneling is also possible between propagating bands. 
However, the main difference is still expected from the 
presence of absence of a flat band (half-integer versus 
integer S), as is also reflected in the different matching 
conditions. 



VIII. POSSIBLE EXPERIMENTS 

In this section, we discuss the expertimental possibil- 
ities to create optical lattices, which would realize the 
spin-S' Dirac equation, and the methods to observe the 
characteristic physical quantities. This section gives a 
brief overview of which concrete protocols have been pro- 
posed in the cold atom literature over the past few years 
- some of them are of course still under active develop- 
ment. Indeed, there is a significant and ongoing experi- 
mental effort devoted to realise lattices with exotic band 
structures, for recent examples on triangle-based lattices, 
see Refsj^i^. Attempts to realise the family of lattices 
proposed here would form part of this endeavour. 

As we have already mentioned, the lattice structure in 
Fig. [T]can be regarded as face centred cubic lattice, which 
in itself - being a Bravail lattice - is relatively straightfor- 
ward to generate. As a first option, this can be created 
by four laser beams at the appropriate angle, realized 
and discussed in detail in Refs<^"— . Although the tri- 
angular layers would be a priori equidistant, a setup like 
that of Ref. already treats the lasers in one [111] direc- 
tion inequivalently from the others, so that the relative 
strengths for intra- and inter-plane hopping need not be 
equal. Making the former much weaker than the latter 
(even without adding further laser beams, relative angles 
and intensities of the beams are tunable) will then yield a 
bandstructure including the presence the flat band, as de- 
scribed above. To achieve the chosen number layers, one 



can e.g. create an optical superlattice in the perpendicu- 
lar direction to the layers^, or by utilizing blue-detuned 
light sheets to terminate the layered structure. Particles 
are then mainly confined to these triangular layers, whose 
number defines 25' -I- 1. 

Second, one can profit from the versatility of a holo- 
graphic mask, enabling arbitrary geometries, to generate 
the desired lattice structure^. Another option is to fol- 
low the steps outlined in Ref<^ and to introduce spin de- 
pendent hopping amplitudes, which in turn also realize 
the desired multiple Dirac-cone structure. 

In terms of observables, the presence of the flat band 
can be revealed by time-of-fiight imaging, since particles 
residing on the flat band remain immobile^, and would 
show up as 'missing' particles. In addition, the number 
of particles (the integral of the density of states) on the 
lattice as a function of the chemical potential could be 
monitored, which a jump around zero energy for integer 
S due to the large degeneracy of the flat band^. The 
particle number per lattice site and physical spin behaves 
close to half fllling as 

N{fi ~ 0) - '^^ ^ ^ = '^■s^'"«^g'^'-^^g^(^) + (45) 

where p(/x) is the DOS, which can be obtained by tak- 
ing the numerical derivative of the experimentally mea- 
sured particle number with respect to fj,. The DOS can 
also conveniently be measured by rf spectroscopy, which 
directly probes the momentum integrated spectral func- 
tion, i.e. the DOS^. The momentum resolved Raman 
spectroscopy can also be used for the same purpose. 

The density-density correlation function, which is 
readily related to the optical conductivity^ can be in- 
vestigated by shot noise measurement, while the optical 
conductivity can directly be probed by the amplitude or 
phase modulation of the optical lattice^. Thanks to the 
modulation, the energy absorption rate or the doublon 
production rate turns out to be directly proportional to 
cr(w). 

To probe the spin-Hall effect, an effective electric held 
should be applied by tilting the lattice along one direc- 
tion, and the detection of the spin current accumula- 
tion through separate imaging of the two different spin 
components^ could reveal the quantization of the spin- 
Hall conductivity, stemming from the underlying topol- 
ogy of the band structure. 

Work in progress on realising artificial gauge fields 
holds the promise to probe the Landau level degeneracy. 
These gauge fields mimic the effect of a real vector poten- 
tial, thus leading to the formation of Landau levels. The 
enhanced degeneracy of the zero energy level should be 
revealed by time-of-flight imaging, as discussed above, or 
by rf spectroscopy. In addition, the Hall current can be 
made visible by driving the system out of equilibrium by 
suddenly changing the trapping potential, and measuring 
the Hall current^. 

Additional methods for detecting topological proper- 
ties, such as quantized Hall conductivity's and chiral 
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edge states^ have been discussed recently in the liter- 
ature and could be generalized to our lattice setup. 

Klein tunneling is expected to be observed in the pres- 
ence of smooth potential barrier, achievable by an ac- 
celerated optical lattice potential^ or by simply tilting 
the lattice^. Both methods would give rise to an ad- 
ditional potential term, varying linearly in one direction 
as V{r) ^ X. A sharp potential barrier is also available 
using the appropriate holographic masks. The charac- 
teristics of Klein tunneling (i.e. perfect transmission at 
given angles) should show up in the measured momentum 
distribution. 

Finally, we mention in passing that photonic crys- 
tals allow the realization of the appropriate lattice 
geometry^, sketched in Fig. [1] The nature of the edge 
states can be probed similarly to Refi^, together with 
the characteristics of Klein tunneling. 

IX. CONCLUSIONS 

In conclusion, we have studied the lattice generaliza- 
tion of the spin- 1/2 Dirac equation of graphene to arbi- 
trary spin. The main difference arises between integer 
and half-integer spins, the former possessing a flat band, 
which is absent for graphene (corresponding to 5 = 1/2). 
As a result, the density of states and the optical conduc- 
tivity are modified, and the topological properties are 
also enriched. We would like to reemphasize the follow- 
ing points. 

First, even in the absence of a perfect S-p Hamiltonian, 
the above multi-cone picture can survive with a different 
ratio of the opening angles between the cones. 

Second, a flat band is expected under general condi- 
tions in the gapless case: any lattice with an odd number 
of layers and with chiral symmetry is expected to have a 
flat band. In the gapped case the flat band is less general, 
but can still be realised in some natural settings, among 
the case of equidistant layers with regular stacking. 

Third, some topological properties depend sensitively 
on the stacking configuration. We have established this 
for the Aharonov-Casher zero-modes in a random mag- 
netic field. Also the Chern numbers of individual bands 
are strongly stacking dependent. 

Fourth, while the wavefunction remains continuous 
across a potential barrier for half-integer pseudospin, 
only its even components remain continuous for integer 
S. 

Last, we would like to emphasize that the lattice struc- 
ture in Fig. [1] simply corresponds to a four-layer slab 
of the face-centred cubic Bravais lattice squeezed in the 
[111] direction. 
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Appendix A: Symmetry conditions for flat bands 

In this appendix we discuss in further detail the differ- 
ent symmetry properties both of the full Hamiltonian in 
Eq. ^ as well as of its low-energy expansion around the 
Dirac points of which the spin S Hamiltonian Eq. ^ 
is a special case. We will show that the existence of the 
fiat band in the different settings can be understood in 
terms of symmetries and that the symmetry conditions 
in the different settings are of different level of generality, 
that is, of different level of robustness to variations in the 
parameters of the model. 

In the gapless case the flat band is protected by a very 
general symmetry — the chiral symmetry — which is made 
possible by the bipartite nature of the lattice, and is 
present under the considered nearest neighbor hopping, 
which respects this bipartitness (bipartiteness will be ex- 
plained below). The flat band is insensitive to arbitrary 
variations of the parameters a^- , for example due to small 
misalignments of the lattice, as this does not change the 
bipartite structure of the Hamiltonian. 

The situation is changed when one introduces any di- 
agonal terms, for example due to intralayer hopping or 
in the form of a gap term like for example A^'^. Such 
terms do not respect the bipartite structure and leads 
to the violation of the the chiral symmetry. However, 
less general symmetries protecting flat bands can still be 
at play. The flat band will no longer be robust to any 
small arbitrary variations in the parameters. However, 
it will still be there for a rather wide class of parameter 
configurations, including fortunately some rather natu- 
ral configurations like equidistant layers. We will now 
discuss these symmetry issues in further detail, starting 
with the chiral symmetry. 

A lattice is bipartite if the sites can be collected into 
two partitions A and B and the Hamiltonian only con- 
tains nonzero matrix elements between the partitions but 
not within each partition (which is the case of the near- 
est neighbor hopping on the honeycomb lattice, which 
gives the free massless Dirac Hamiltonian). In our mul- 
tilayer generalization the partitions are the odd and the 
even layers, respectively. Chiral symmetry is present as 
the only terms in the Hamiltonian are the hoppings be- 
tween adjacent layers. In contrast, hopping within the 
layers of triangular lattices would be one source of diago- 
nal terms in the Hamiltonian Eq. 3] A transverse electric 
field would be another source. Such diagonal terms break 
the chiral symmetry. 
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We are now going to discuss on a more mathematical level the conjugation symmetry properties of the model 

Hamiltonian 



Hs = t 



Ao aoi/(k) 
«Si/*(k) Ai 
alJ*{k)^ 






ai2/(k)^ 

a2S-i,25/(k) 

45-l,25/*(k) As 



(Al) 



When the gap terms are zero A.; = 0, the Hamiltonian 
possesses a chiral symmetry (CS): I]H{'k)T,^ = —H(k), 
where E is a unitary matrix. The CS conjugates the 
spectrum. In our case we have E = diag(l, —1, 1, —1, . . .) 
generalizing E = cr^ that conjugates the gapless S ~ ^ 
Dirac spectrum. In the case of an odd number of bands 
the CS has to map the middlemost band onto itself. 
Furthermore, since the CS conjugates the spectrum at 
each k separately, each (crystal) momentum eigenstate 
in the middle band has to map onto itself. Chiral sym- 
metry guarantees thus a zero-energy state at each k — a 
flat band. This should be contrasted with particle-hole 
symmetry (PHS) UpKH*{k)Ul^ = -H{-k), also pos- 
sessed by the same Hamiltonian. PHS conjugates the 
band structure, however, does not do so at each k inde- 
pendently wherefore PHS alone is not enough to guaran- 
tee a flat bandit 

The chiral symmetry is not only independent of the 
precise values of the parameters a^j's, but also of the 
stacking. As the chiral symmetry only relies on bipar- 
titeness, it is also independent to changes in the stack- 
ing with the consequent individual flipping of chiralities 
in the tunneling between the involved layers, as for the 
Hamiltonian Eq. ([T0| . 

For the gapped case the chiral symmetry is lost and 
there are no general symmetries to guarantee a flat 
band. Nonetheless, the following conditions are still suf- 
ficient for the appearance of a flat band in the gen- 
eral case with a nonzero gap. The diagonal gap term 
diag(Ao, Al, . . . , A2s) must satisfy the anti-symmetric 
property A25 = — Aq, A2S-1 = — Ai etc. At the 
same time, the off-diagonal elements with the a param- 
eters must satisfy a symmetric property in their mod- 



uH: |a2S-i,s| = |aoi|, |a2S-2,s-i| = |ai2| etc., while 
the phases can be chosen arbitrarily. For this reason, 
the flipping of chiralities between layers does not change 
anything since |/*| = |/|. 

Even though it is difficult to interpret these general 
conditions in terms of symmetries, these can be con- 
nected, in some special cases, to a symmetry property 
specific to the Dirac cone, that is, they are present only 
in the low-energy approximation for the system close to 
half-filling. For example, if one imposes the stronger con- 
dition q;2S-i,S = aoi, "25-2,5-1 = ai2, •■• and a sym- 
metric configuration of flipped chiralities, then the fiat 
bands can be seen as a consequence of the emergent con- 
jugation property YHg{p)Y^ = —Hs{p), with p being 
the small momentum with respect to a Dirac point and 
with Y generalizing (7 y dbS Oil matrix with the alternat- 
ing pattern (— i, i, — i, i, . . .) on the anti-diagonal of the 
matrix (i.e., the diagonal joining bottom left and top 
right entries). It formally looks like a PH symmetry com- 
bined with inversion symmetry H{—p) = H{p), however, 
with momenta not inverted around k = as the phys- 
ical PH symmetry requires but around the Dirac point 
p = k - k+ = 0. 

As mentioned in the beginning of the section, some 
of the specific parameter configurations are rather nat- 
ural. For equidistant layers, the gap parameter is ASz 
and can be created by the electrostatic potential of a 
uniform transverse electric field. The a's would all have 
the same modulus for equidistant layers and therefore 
also fulfill the above conditions, although they would not 
correspond to the values pertinent to Sx and Sy. Thus, 
we conclude that equidistant layers for certain symmetric 
stackings are enough to guarantee flat bands also in the 
presence of a gap term generated by a transverse uniform 
electric field. 
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